<!DOCTYPE html>
<HTML lang = "en">
<HEAD>
  <meta charset="UTF-8"/>
  <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
  <title>The Different Flavors of Parallelism</title>
  

  <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]},
      TeX: { equationNumbers: { autoNumber: "AMS" } }
    });
  </script>

  <script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
  </script>

  
<style>
pre.hljl {
    border: 1px solid #ccc;
    margin: 5px;
    padding: 5px;
    overflow-x: auto;
    color: rgb(68,68,68); background-color: rgb(251,251,251); }
pre.hljl > span.hljl-t { }
pre.hljl > span.hljl-w { }
pre.hljl > span.hljl-e { }
pre.hljl > span.hljl-eB { }
pre.hljl > span.hljl-o { }
pre.hljl > span.hljl-k { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kc { color: rgb(59,151,46); font-style: italic; }
pre.hljl > span.hljl-kd { color: rgb(214,102,97); font-style: italic; }
pre.hljl > span.hljl-kn { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kp { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kr { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kt { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-n { }
pre.hljl > span.hljl-na { }
pre.hljl > span.hljl-nb { }
pre.hljl > span.hljl-nbp { }
pre.hljl > span.hljl-nc { }
pre.hljl > span.hljl-ncB { }
pre.hljl > span.hljl-nd { color: rgb(214,102,97); }
pre.hljl > span.hljl-ne { }
pre.hljl > span.hljl-neB { }
pre.hljl > span.hljl-nf { color: rgb(66,102,213); }
pre.hljl > span.hljl-nfm { color: rgb(66,102,213); }
pre.hljl > span.hljl-np { }
pre.hljl > span.hljl-nl { }
pre.hljl > span.hljl-nn { }
pre.hljl > span.hljl-no { }
pre.hljl > span.hljl-nt { }
pre.hljl > span.hljl-nv { }
pre.hljl > span.hljl-nvc { }
pre.hljl > span.hljl-nvg { }
pre.hljl > span.hljl-nvi { }
pre.hljl > span.hljl-nvm { }
pre.hljl > span.hljl-l { }
pre.hljl > span.hljl-ld { color: rgb(148,91,176); font-style: italic; }
pre.hljl > span.hljl-s { color: rgb(201,61,57); }
pre.hljl > span.hljl-sa { color: rgb(201,61,57); }
pre.hljl > span.hljl-sb { color: rgb(201,61,57); }
pre.hljl > span.hljl-sc { color: rgb(201,61,57); }
pre.hljl > span.hljl-sd { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdB { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdC { color: rgb(201,61,57); }
pre.hljl > span.hljl-se { color: rgb(59,151,46); }
pre.hljl > span.hljl-sh { color: rgb(201,61,57); }
pre.hljl > span.hljl-si { }
pre.hljl > span.hljl-so { color: rgb(201,61,57); }
pre.hljl > span.hljl-sr { color: rgb(201,61,57); }
pre.hljl > span.hljl-ss { color: rgb(201,61,57); }
pre.hljl > span.hljl-ssB { color: rgb(201,61,57); }
pre.hljl > span.hljl-nB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nbB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nfB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nh { color: rgb(59,151,46); }
pre.hljl > span.hljl-ni { color: rgb(59,151,46); }
pre.hljl > span.hljl-nil { color: rgb(59,151,46); }
pre.hljl > span.hljl-noB { color: rgb(59,151,46); }
pre.hljl > span.hljl-oB { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-ow { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-p { }
pre.hljl > span.hljl-c { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-ch { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cm { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cp { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cpB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cs { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-csB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-g { }
pre.hljl > span.hljl-gd { }
pre.hljl > span.hljl-ge { }
pre.hljl > span.hljl-geB { }
pre.hljl > span.hljl-gh { }
pre.hljl > span.hljl-gi { }
pre.hljl > span.hljl-go { }
pre.hljl > span.hljl-gp { }
pre.hljl > span.hljl-gs { }
pre.hljl > span.hljl-gsB { }
pre.hljl > span.hljl-gt { }
</style>



  <style type="text/css">
  @font-face {
  font-style: normal;
  font-weight: 300;
}
@font-face {
  font-style: normal;
  font-weight: 400;
}
@font-face {
  font-style: normal;
  font-weight: 600;
}
html {
  font-family: sans-serif; /* 1 */
  -ms-text-size-adjust: 100%; /* 2 */
  -webkit-text-size-adjust: 100%; /* 2 */
}
body {
  margin: 0;
}
article,
aside,
details,
figcaption,
figure,
footer,
header,
hgroup,
main,
menu,
nav,
section,
summary {
  display: block;
}
audio,
canvas,
progress,
video {
  display: inline-block; /* 1 */
  vertical-align: baseline; /* 2 */
}
audio:not([controls]) {
  display: none;
  height: 0;
}
[hidden],
template {
  display: none;
}
a:active,
a:hover {
  outline: 0;
}
abbr[title] {
  border-bottom: 1px dotted;
}
b,
strong {
  font-weight: bold;
}
dfn {
  font-style: italic;
}
h1 {
  font-size: 2em;
  margin: 0.67em 0;
}
mark {
  background: #ff0;
  color: #000;
}
small {
  font-size: 80%;
}
sub,
sup {
  font-size: 75%;
  line-height: 0;
  position: relative;
  vertical-align: baseline;
}
sup {
  top: -0.5em;
}
sub {
  bottom: -0.25em;
}
img {
  border: 0;
}
svg:not(:root) {
  overflow: hidden;
}
figure {
  margin: 1em 40px;
}
hr {
  -moz-box-sizing: content-box;
  box-sizing: content-box;
  height: 0;
}
pre {
  overflow: auto;
}
code,
kbd,
pre,
samp {
  font-family: monospace, monospace;
  font-size: 1em;
}
button,
input,
optgroup,
select,
textarea {
  color: inherit; /* 1 */
  font: inherit; /* 2 */
  margin: 0; /* 3 */
}
button {
  overflow: visible;
}
button,
select {
  text-transform: none;
}
button,
html input[type="button"], /* 1 */
input[type="reset"],
input[type="submit"] {
  -webkit-appearance: button; /* 2 */
  cursor: pointer; /* 3 */
}
button[disabled],
html input[disabled] {
  cursor: default;
}
button::-moz-focus-inner,
input::-moz-focus-inner {
  border: 0;
  padding: 0;
}
input {
  line-height: normal;
}
input[type="checkbox"],
input[type="radio"] {
  box-sizing: border-box; /* 1 */
  padding: 0; /* 2 */
}
input[type="number"]::-webkit-inner-spin-button,
input[type="number"]::-webkit-outer-spin-button {
  height: auto;
}
input[type="search"] {
  -webkit-appearance: textfield; /* 1 */
  -moz-box-sizing: content-box;
  -webkit-box-sizing: content-box; /* 2 */
  box-sizing: content-box;
}
input[type="search"]::-webkit-search-cancel-button,
input[type="search"]::-webkit-search-decoration {
  -webkit-appearance: none;
}
fieldset {
  border: 1px solid #c0c0c0;
  margin: 0 2px;
  padding: 0.35em 0.625em 0.75em;
}
legend {
  border: 0; /* 1 */
  padding: 0; /* 2 */
}
textarea {
  overflow: auto;
}
optgroup {
  font-weight: bold;
}
table {
  font-family: monospace, monospace;
  font-size : 0.8em;
  border-collapse: collapse;
  border-spacing: 0;
}
td,
th {
  padding: 0;
}
thead th {
    border-bottom: 1px solid black;
    background-color: white;
}
tr:nth-child(odd){
  background-color: rgb(248,248,248);
}


/*
* Skeleton V2.0.4
* Copyright 2014, Dave Gamache
* www.getskeleton.com
* Free to use under the MIT license.
* http://www.opensource.org/licenses/mit-license.php
* 12/29/2014
*/
.container {
  position: relative;
  width: 100%;
  max-width: 960px;
  margin: 0 auto;
  padding: 0 20px;
  box-sizing: border-box; }
.column,
.columns {
  width: 100%;
  float: left;
  box-sizing: border-box; }
@media (min-width: 400px) {
  .container {
    width: 85%;
    padding: 0; }
}
@media (min-width: 550px) {
  .container {
    width: 80%; }
  .column,
  .columns {
    margin-left: 4%; }
  .column:first-child,
  .columns:first-child {
    margin-left: 0; }

  .one.column,
  .one.columns                    { width: 4.66666666667%; }
  .two.columns                    { width: 13.3333333333%; }
  .three.columns                  { width: 22%;            }
  .four.columns                   { width: 30.6666666667%; }
  .five.columns                   { width: 39.3333333333%; }
  .six.columns                    { width: 48%;            }
  .seven.columns                  { width: 56.6666666667%; }
  .eight.columns                  { width: 65.3333333333%; }
  .nine.columns                   { width: 74.0%;          }
  .ten.columns                    { width: 82.6666666667%; }
  .eleven.columns                 { width: 91.3333333333%; }
  .twelve.columns                 { width: 100%; margin-left: 0; }

  .one-third.column               { width: 30.6666666667%; }
  .two-thirds.column              { width: 65.3333333333%; }

  .one-half.column                { width: 48%; }

  /* Offsets */
  .offset-by-one.column,
  .offset-by-one.columns          { margin-left: 8.66666666667%; }
  .offset-by-two.column,
  .offset-by-two.columns          { margin-left: 17.3333333333%; }
  .offset-by-three.column,
  .offset-by-three.columns        { margin-left: 26%;            }
  .offset-by-four.column,
  .offset-by-four.columns         { margin-left: 34.6666666667%; }
  .offset-by-five.column,
  .offset-by-five.columns         { margin-left: 43.3333333333%; }
  .offset-by-six.column,
  .offset-by-six.columns          { margin-left: 52%;            }
  .offset-by-seven.column,
  .offset-by-seven.columns        { margin-left: 60.6666666667%; }
  .offset-by-eight.column,
  .offset-by-eight.columns        { margin-left: 69.3333333333%; }
  .offset-by-nine.column,
  .offset-by-nine.columns         { margin-left: 78.0%;          }
  .offset-by-ten.column,
  .offset-by-ten.columns          { margin-left: 86.6666666667%; }
  .offset-by-eleven.column,
  .offset-by-eleven.columns       { margin-left: 95.3333333333%; }

  .offset-by-one-third.column,
  .offset-by-one-third.columns    { margin-left: 34.6666666667%; }
  .offset-by-two-thirds.column,
  .offset-by-two-thirds.columns   { margin-left: 69.3333333333%; }

  .offset-by-one-half.column,
  .offset-by-one-half.columns     { margin-left: 52%; }

}
html {
  font-size: 62.5%; }
body {
  font-size: 1.5em; /* currently ems cause chrome bug misinterpreting rems on body element */
  line-height: 1.6;
  font-weight: 400;
  font-family: "Raleway", "HelveticaNeue", "Helvetica Neue", Helvetica, Arial, sans-serif;
  color: #222; }
h1, h2, h3, h4, h5, h6 {
  margin-top: 0;
  margin-bottom: 2rem;
  font-weight: 300; }
h1 { font-size: 3.6rem; line-height: 1.2;  letter-spacing: -.1rem;}
h2 { font-size: 3.4rem; line-height: 1.25; letter-spacing: -.1rem; }
h3 { font-size: 3.2rem; line-height: 1.3;  letter-spacing: -.1rem; }
h4 { font-size: 2.8rem; line-height: 1.35; letter-spacing: -.08rem; }
h5 { font-size: 2.4rem; line-height: 1.5;  letter-spacing: -.05rem; }
h6 { font-size: 1.5rem; line-height: 1.6;  letter-spacing: 0; }

p {
  margin-top: 0; }
a {
  color: #1EAEDB; }
a:hover {
  color: #0FA0CE; }
.button,
button,
input[type="submit"],
input[type="reset"],
input[type="button"] {
  display: inline-block;
  height: 38px;
  padding: 0 30px;
  color: #555;
  text-align: center;
  font-size: 11px;
  font-weight: 600;
  line-height: 38px;
  letter-spacing: .1rem;
  text-transform: uppercase;
  text-decoration: none;
  white-space: nowrap;
  background-color: transparent;
  border-radius: 4px;
  border: 1px solid #bbb;
  cursor: pointer;
  box-sizing: border-box; }
.button:hover,
button:hover,
input[type="submit"]:hover,
input[type="reset"]:hover,
input[type="button"]:hover,
.button:focus,
button:focus,
input[type="submit"]:focus,
input[type="reset"]:focus,
input[type="button"]:focus {
  color: #333;
  border-color: #888;
  outline: 0; }
.button.button-primary,
button.button-primary,
input[type="submit"].button-primary,
input[type="reset"].button-primary,
input[type="button"].button-primary {
  color: #FFF;
  background-color: #33C3F0;
  border-color: #33C3F0; }
.button.button-primary:hover,
button.button-primary:hover,
input[type="submit"].button-primary:hover,
input[type="reset"].button-primary:hover,
input[type="button"].button-primary:hover,
.button.button-primary:focus,
button.button-primary:focus,
input[type="submit"].button-primary:focus,
input[type="reset"].button-primary:focus,
input[type="button"].button-primary:focus {
  color: #FFF;
  background-color: #1EAEDB;
  border-color: #1EAEDB; }
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea,
select {
  height: 38px;
  padding: 6px 10px; /* The 6px vertically centers text on FF, ignored by Webkit */
  background-color: #fff;
  border: 1px solid #D1D1D1;
  border-radius: 4px;
  box-shadow: none;
  box-sizing: border-box; }
/* Removes awkward default styles on some inputs for iOS */
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea {
  -webkit-appearance: none;
     -moz-appearance: none;
          appearance: none; }
textarea {
  min-height: 65px;
  padding-top: 6px;
  padding-bottom: 6px; }
input[type="email"]:focus,
input[type="number"]:focus,
input[type="search"]:focus,
input[type="text"]:focus,
input[type="tel"]:focus,
input[type="url"]:focus,
input[type="password"]:focus,
textarea:focus,
select:focus {
  border: 1px solid #33C3F0;
  outline: 0; }
label,
legend {
  display: block;
  margin-bottom: .5rem;
  font-weight: 600; }
fieldset {
  padding: 0;
  border-width: 0; }
input[type="checkbox"],
input[type="radio"] {
  display: inline; }
label > .label-body {
  display: inline-block;
  margin-left: .5rem;
  font-weight: normal; }
ul {
  list-style: circle; }
ol {
  list-style: decimal; }
ul ul,
ul ol,
ol ol,
ol ul {
  margin: 1.5rem 0 1.5rem 3rem;
  font-size: 90%; }
li > p {margin : 0;}
th,
td {
  padding: 12px 15px;
  text-align: left;
  border-bottom: 1px solid #E1E1E1; }
th:first-child,
td:first-child {
  padding-left: 0; }
th:last-child,
td:last-child {
  padding-right: 0; }
button,
.button {
  margin-bottom: 1rem; }
input,
textarea,
select,
fieldset {
  margin-bottom: 1.5rem; }
pre,
blockquote,
dl,
figure,
table,
p,
ul,
ol,
form {
  margin-bottom: 1.0rem; }
.u-full-width {
  width: 100%;
  box-sizing: border-box; }
.u-max-full-width {
  max-width: 100%;
  box-sizing: border-box; }
.u-pull-right {
  float: right; }
.u-pull-left {
  float: left; }
hr {
  margin-top: 3rem;
  margin-bottom: 3.5rem;
  border-width: 0;
  border-top: 1px solid #E1E1E1; }
.container:after,
.row:after,
.u-cf {
  content: "";
  display: table;
  clear: both; }

pre {
  display: block;
  padding: 9.5px;
  margin: 0 0 10px;
  font-size: 13px;
  line-height: 1.42857143;
  word-break: break-all;
  word-wrap: break-word;
  border: 1px solid #ccc;
  border-radius: 4px;
}

pre.hljl {
  margin: 0 0 10px;
  display: block;
  background: #f5f5f5;
  border-radius: 4px;
  padding : 5px;
}

pre.output {
  background: #ffffff;
}

pre.code {
  background: #ffffff;
}

pre.julia-error {
  color : red
}

code,
kbd,
pre,
samp {
  font-family: Menlo, Monaco, Consolas, "Courier New", monospace;
  font-size: 0.9em;
}


@media (min-width: 400px) {}
@media (min-width: 550px) {}
@media (min-width: 750px) {}
@media (min-width: 1000px) {}
@media (min-width: 1200px) {}

h1.title {margin-top : 20px}
img {max-width : 100%}
div.title {text-align: center;}

  </style>
</HEAD>

<BODY>
  <div class ="container">
    <div class = "row">
      <div class = "col-md-12 twelve columns">
        <div class="title">
          <h1 class="title">The Different Flavors of Parallelism</h1>
          <h5>Chris Rackauckas</h5>
          <h5>September 25th, 2020</h5>
        </div>

        <p>Now that you are aware of the basics of parallel computing, let&#39;s give a high level overview of the differences between different modes of parallelism.</p>
<h2>Lowest Level: SIMD</h2>
<p>Recall SIMD, the idea that processors can run multiple commands simulataniously on specially structured data. &quot;Single Instruction Multiple Data&quot;. SIMD is parallelism within a single core.</p>
<h3>High Level Idea of SIMD</h3>
<p>Calculations can occur in parallel in the processor if there is sufficient structure in the computation.</p>
<h3>How to do SIMD</h3>
<p>The simplest way to do SIMD is simply to make sure that your values are aligned. If they are, then great, LLVM&#39;s autovectorizer pass has a good chance of automatic vectorization &#40;in the world of computing, &quot;SIMD&quot; is synonymous with vectorization since it is taking specific values and instead computing on small vectors. That is not to be confused with &quot;vectorization&quot; in the sense of Python/R/MATLAB, which is a programming style which prefers using C-defined primitive functions, like broadcast or matrix multiplication&#41;.</p>
<p>You can check for auto-vectorization inside of the LLVM IR by looking for statements like:</p>
<pre><code>&#37;wide.load24 &#61; load &lt;4 x double&gt;, &lt;4 x double&gt; addrspac&#40;13&#41;* &#37;46, align 8
; └
; ┌ @ float.jl:395 within &#96;&#43;&#39;
&#37;47 &#61; fadd &lt;4 x double&gt; &#37;wide.load, &#37;wide.load24</code></pre>
<p>which means that 4 additions are happening simultaniously. The amount of vectorization is heavily dependent on your architecture. The ancient form of SIMD, the SSE&#40;2&#41; instructions, required that your data was aligned. Now there&#39;s a bit more leeway, but generally it holds that making your the data you&#39;re trying to SIMD over is aligned. Thus there can be major differences in computing using a <em>struct of array</em> format instead of an <em>arrays of structs</em> format. For example:</p>


<pre class='hljl'>
<span class='hljl-k'>struct</span><span class='hljl-t'> </span><span class='hljl-n'>MyComplex</span><span class='hljl-t'>
  </span><span class='hljl-n'>real</span><span class='hljl-oB'>::</span><span class='hljl-n'>Float64</span><span class='hljl-t'>
  </span><span class='hljl-n'>imag</span><span class='hljl-oB'>::</span><span class='hljl-n'>Float64</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>arr</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nf'>MyComplex</span><span class='hljl-p'>(</span><span class='hljl-nf'>rand</span><span class='hljl-p'>(),</span><span class='hljl-nf'>rand</span><span class='hljl-p'>())</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>100</span><span class='hljl-p'>]</span>
</pre>


<pre class="output">
100-element Array&#123;Main.##WeaveSandBox#406.MyComplex,1&#125;:
 Main.##WeaveSandBox#406.MyComplex&#40;0.8097526244493014, 0.7421937575992628&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.13110792466742582, 0.15960559660680018
&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.5342419774726372, 0.3338732796442523&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.6376540382212976, 0.1569913652842596&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.12060395270536128, 0.9419349640236092&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.9407425560872236, 0.7549489482002789&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.2507206100986332, 0.2408940656163221&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.7089307125188724, 0.777669384190381&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.6534826600645793, 0.6499210868969258&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.5297240732394792, 0.7910214341072723&#41;
 ⋮
 Main.##WeaveSandBox#406.MyComplex&#40;0.014577544502819562, 0.7485575815519052
&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.2957281165123682, 0.32356309079613155&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.9998411644215583, 0.8591886629678238&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.6025356311401455, 0.5521040680074492&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.4185874998482655, 0.41723979114768106&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.6343517522572546, 0.20332920724512582&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.4187947523552482, 0.5534240872727254&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.3736942345513681, 0.5524908964899642&#41;
 Main.##WeaveSandBox#406.MyComplex&#40;0.25105642946769047, 0.6786745580088376&#41;
</pre>


<p>is represented in memory as</p>
<pre><code>&#91;real1,imag1,real2,imag2,...&#93;</code></pre>
<p>while the struct of array formats are</p>


<pre class='hljl'>
<span class='hljl-k'>struct</span><span class='hljl-t'> </span><span class='hljl-n'>MyComplexes</span><span class='hljl-t'>
  </span><span class='hljl-n'>real</span><span class='hljl-oB'>::</span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-n'>Float64</span><span class='hljl-p'>}</span><span class='hljl-t'>
  </span><span class='hljl-n'>imag</span><span class='hljl-oB'>::</span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-n'>Float64</span><span class='hljl-p'>}</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>arr2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>MyComplexes</span><span class='hljl-p'>(</span><span class='hljl-nf'>rand</span><span class='hljl-p'>(</span><span class='hljl-ni'>100</span><span class='hljl-p'>),</span><span class='hljl-nf'>rand</span><span class='hljl-p'>(</span><span class='hljl-ni'>100</span><span class='hljl-p'>))</span>
</pre>


<pre class="output">
Main.##WeaveSandBox#406.MyComplexes&#40;&#91;0.22551923261187778, 0.026566259011969
74, 0.8113637753837821, 0.5255795986450098, 0.7620987825532188, 0.444655447
0643247, 0.4315190560085518, 0.35440840549215813, 0.0029381746761620864, 0.
3807316833053711  …  0.2747317708028554, 0.7393153412285369, 0.256618072905
38907, 0.2934818668478125, 0.8140976859891642, 0.6330027742505249, 0.232768
2560032629, 0.561810292923401, 0.02584282653153447, 0.8183914022941403&#93;, &#91;0
.754509086722589, 0.062597821009311, 0.9173163667357256, 0.8054991166605592
, 0.5117437953535597, 0.6134087820806926, 0.9926351239833742, 0.46494820009
49148, 0.5218766840397984, 0.4201587652158061  …  0.8970547855806963, 0.996
4322566447046, 0.1681551551274456, 0.911100131963601, 0.4979583988274223, 0
.19617500522829112, 0.4563113262874552, 0.6146552580989488, 0.8835704379028
788, 0.7641379717749286&#93;&#41;
</pre>


<p>Now let&#39;s check what happens when we perform a reduction:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>InteractiveUtils</span><span class='hljl-t'>
</span><span class='hljl-n'>Base</span><span class='hljl-oB'>.:+</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-oB'>::</span><span class='hljl-n'>MyComplex</span><span class='hljl-p'>,</span><span class='hljl-n'>y</span><span class='hljl-oB'>::</span><span class='hljl-n'>MyComplex</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>MyComplex</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-oB'>.</span><span class='hljl-n'>real</span><span class='hljl-oB'>+</span><span class='hljl-n'>y</span><span class='hljl-oB'>.</span><span class='hljl-n'>real</span><span class='hljl-p'>,</span><span class='hljl-n'>x</span><span class='hljl-oB'>.</span><span class='hljl-n'>imag</span><span class='hljl-oB'>+</span><span class='hljl-n'>y</span><span class='hljl-oB'>.</span><span class='hljl-n'>imag</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>Base</span><span class='hljl-oB'>.:/</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-oB'>::</span><span class='hljl-n'>MyComplex</span><span class='hljl-p'>,</span><span class='hljl-n'>y</span><span class='hljl-oB'>::</span><span class='hljl-n'>Int</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>MyComplex</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-oB'>.</span><span class='hljl-n'>real</span><span class='hljl-oB'>/</span><span class='hljl-n'>y</span><span class='hljl-p'>,</span><span class='hljl-n'>x</span><span class='hljl-oB'>.</span><span class='hljl-n'>imag</span><span class='hljl-oB'>/</span><span class='hljl-n'>y</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>average</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-oB'>::</span><span class='hljl-nf'>Vector</span><span class='hljl-p'>{</span><span class='hljl-n'>MyComplex</span><span class='hljl-p'>})</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>sum</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-p'>)</span><span class='hljl-oB'>/</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nd'>@code_llvm</span><span class='hljl-t'> </span><span class='hljl-nf'>average</span><span class='hljl-p'>(</span><span class='hljl-n'>arr</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
;  @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.jm
d:5 within &#96;average&#39;
; Function Attrs: uwtable
define void @julia_average_11839&#40;&#91;2 x double&#93;* noalias nocapture sret, &#37;jl_
value_t* nonnull align 16 dereferenceable&#40;40&#41;&#41; #0 &#123;
top:
  &#37;2 &#61; alloca &#37;jl_value_t*, i32 4
  &#37;3 &#61; alloca &#91;2 x double&#93;, align 8
; ┌ @ reducedim.jl:719 within &#96;sum&#39;
; │┌ @ reducedim.jl:719 within &#96;#sum#627&#39;
; ││┌ @ reducedim.jl:723 within &#96;_sum&#39; @ reducedim.jl:724
; │││┌ @ reducedim.jl:310 within &#96;mapreduce&#39;
; ││││┌ @ reducedim.jl:310 within &#96;#mapreduce#620&#39;
; │││││┌ @ reducedim.jl:318 within &#96;_mapreduce_dim&#39;
; ││││││┌ @ reduce.jl:397 within &#96;_mapreduce&#39;
; │││││││┌ @ indices.jl:453 within &#96;LinearIndices&#39;
; ││││││││┌ @ abstractarray.jl:75 within &#96;axes&#39;
; │││││││││┌ @ array.jl:155 within &#96;size&#39;
            &#37;4 &#61; bitcast &#37;jl_value_t* &#37;1 to &#37;jl_value_t**
            &#37;5 &#61; getelementptr inbounds &#37;jl_value_t*, &#37;jl_value_t** &#37;4, i64
 3
            &#37;6 &#61; bitcast &#37;jl_value_t** &#37;5 to i64*
            &#37;7 &#61; load i64, i64* &#37;6, align 8
; │││││││││└
; │││││││││┌ @ tuple.jl:157 within &#96;map&#39;
; ││││││││││┌ @ range.jl:326 within &#96;OneTo&#39; @ range.jl:317
; │││││││││││┌ @ promotion.jl:409 within &#96;max&#39;
              &#37;8 &#61; icmp sgt i64 &#37;7, 0
              &#37;9 &#61; select i1 &#37;8, i64 &#37;7, i64 0
; │││││││└└└└└
; │││││││ @ reduce.jl:399 within &#96;_mapreduce&#39;
         br i1 &#37;8, label &#37;L11, label &#37;L9

L9:                                               ; preds &#61; &#37;top
; │││││││ @ reduce.jl:400 within &#96;_mapreduce&#39;
         &#37;10 &#61; getelementptr &#37;jl_value_t*, &#37;jl_value_t** &#37;2, i32 0
         store &#37;jl_value_t* inttoptr &#40;i64 191572544 to &#37;jl_value_t*&#41;, &#37;jl_v
alue_t** &#37;10
         &#37;11 &#61; getelementptr &#37;jl_value_t*, &#37;jl_value_t** &#37;2, i32 1
         store &#37;jl_value_t* inttoptr &#40;i64 218766144 to &#37;jl_value_t*&#41;, &#37;jl_v
alue_t** &#37;11
         &#37;12 &#61; getelementptr &#37;jl_value_t*, &#37;jl_value_t** &#37;2, i32 2
         store &#37;jl_value_t* &#37;1, &#37;jl_value_t** &#37;12
         &#37;13 &#61; getelementptr &#37;jl_value_t*, &#37;jl_value_t** &#37;2, i32 3
         store &#37;jl_value_t* inttoptr &#40;i64 183661808 to &#37;jl_value_t*&#41;, &#37;jl_v
alue_t** &#37;13
         &#37;14 &#61; call nonnull &#37;jl_value_t* @jl_invoke&#40;&#37;jl_value_t* inttoptr &#40;
i64 191682656 to &#37;jl_value_t*&#41;, &#37;jl_value_t** &#37;2, i32 4, &#37;jl_value_t* intto
ptr &#40;i64 954352528 to &#37;jl_value_t*&#41;&#41;
         call void @llvm.trap&#40;&#41;
         unreachable

L11:                                              ; preds &#61; &#37;top
; │││││││ @ reduce.jl:401 within &#96;_mapreduce&#39;
; │││││││┌ @ promotion.jl:398 within &#96;&#61;&#61;&#39;
          &#37;15 &#61; icmp eq i64 &#37;9, 1
; │││││││└
         br i1 &#37;15, label &#37;L13, label &#37;L15

L13:                                              ; preds &#61; &#37;L11
; │││││││ @ reduce.jl:402 within &#96;_mapreduce&#39;
; │││││││┌ @ array.jl:809 within &#96;getindex&#39;
          &#37;16 &#61; bitcast &#37;jl_value_t* &#37;1 to &#91;2 x double&#93;**
          &#37;17 &#61; load &#91;2 x double&#93;*, &#91;2 x double&#93;** &#37;16, align 8
          &#37;.elt &#61; getelementptr inbounds &#91;2 x double&#93;, &#91;2 x double&#93;* &#37;17, i
64 0, i64 0
          &#37;18 &#61; bitcast double* &#37;.elt to &lt;2 x double&gt;*
          &#37;19 &#61; load &lt;2 x double&gt;, &lt;2 x double&gt;* &#37;18, align 8
; │││││││└
; │││││││ @ reduce.jl:403 within &#96;_mapreduce&#39;
         br label &#37;L47

L15:                                              ; preds &#61; &#37;L11
; │││││││ @ reduce.jl:404 within &#96;_mapreduce&#39;
; │││││││┌ @ int.jl:82 within &#96;&lt;&#39;
          &#37;20 &#61; icmp sgt i64 &#37;9, 15
; │││││││└
         br i1 &#37;20, label &#37;L43, label &#37;L17

L17:                                              ; preds &#61; &#37;L15
; │││││││ @ reduce.jl:406 within &#96;_mapreduce&#39;
; │││││││┌ @ array.jl:809 within &#96;getindex&#39;
          &#37;21 &#61; bitcast &#37;jl_value_t* &#37;1 to &#91;2 x double&#93;**
          &#37;22 &#61; load &#91;2 x double&#93;*, &#91;2 x double&#93;** &#37;21, align 8
          &#37;.elt60 &#61; getelementptr inbounds &#91;2 x double&#93;, &#91;2 x double&#93;* &#37;22,
 i64 0, i64 0
          &#37;23 &#61; bitcast double* &#37;.elt60 to &lt;2 x double&gt;*
          &#37;24 &#61; load &lt;2 x double&gt;, &lt;2 x double&gt;* &#37;23, align 8
; │││││││└
; │││││││ @ reduce.jl:407 within &#96;_mapreduce&#39;
; │││││││┌ @ array.jl:809 within &#96;getindex&#39;
          &#37;.elt64 &#61; getelementptr inbounds &#91;2 x double&#93;, &#91;2 x double&#93;* &#37;22,
 i64 1, i64 0
          &#37;25 &#61; bitcast double* &#37;.elt64 to &lt;2 x double&gt;*
          &#37;26 &#61; load &lt;2 x double&gt;, &lt;2 x double&gt;* &#37;25, align 8
; │││││││└
; │││││││ @ reduce.jl:408 within &#96;_mapreduce&#39;
; │││││││┌ @ reduce.jl:24 within &#96;add_sum&#39;
; ││││││││┌ @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_paral
lelism.jmd:3 within &#96;&#43;&#39; @ float.jl:401
           &#37;27 &#61; fadd &lt;2 x double&gt; &#37;24, &#37;26
; │││││││└└
; │││││││ @ reduce.jl:409 within &#96;_mapreduce&#39;
; │││││││┌ @ int.jl:82 within &#96;&lt;&#39;
          &#37;28 &#61; icmp sgt i64 &#37;9, 2
; │││││││└
         br i1 &#37;28, label &#37;L34, label &#37;L47

L34:                                              ; preds &#61; &#37;L34, &#37;L17
         &#37;value_phi74 &#61; phi i64 &#91; &#37;30, &#37;L34 &#93;, &#91; 2, &#37;L17 &#93;
         &#37;29 &#61; phi &lt;2 x double&gt; &#91; &#37;33, &#37;L34 &#93;, &#91; &#37;27, &#37;L17 &#93;
; │││││││ @ reduce.jl:410 within &#96;_mapreduce&#39;
; │││││││┌ @ int.jl:86 within &#96;&#43;&#39;
          &#37;30 &#61; add nuw nsw i64 &#37;value_phi74, 1
; │││││││└
; │││││││┌ @ array.jl:809 within &#96;getindex&#39;
          &#37;.elt68 &#61; getelementptr inbounds &#91;2 x double&#93;, &#91;2 x double&#93;* &#37;22,
 i64 &#37;value_phi74, i64 0
          &#37;31 &#61; bitcast double* &#37;.elt68 to &lt;2 x double&gt;*
          &#37;32 &#61; load &lt;2 x double&gt;, &lt;2 x double&gt;* &#37;31, align 8
; │││││││└
; │││││││ @ reduce.jl:411 within &#96;_mapreduce&#39;
; │││││││┌ @ reduce.jl:24 within &#96;add_sum&#39;
; ││││││││┌ @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_paral
lelism.jmd:3 within &#96;&#43;&#39; @ float.jl:401
           &#37;33 &#61; fadd &lt;2 x double&gt; &#37;29, &#37;32
; │││││││└└
; │││││││ @ reduce.jl:409 within &#96;_mapreduce&#39;
; │││││││┌ @ int.jl:82 within &#96;&lt;&#39;
          &#37;34 &#61; icmp ult i64 &#37;30, &#37;9
; │││││││└
         br i1 &#37;34, label &#37;L34, label &#37;L47

L43:                                              ; preds &#61; &#37;L15
; │││││││ @ reduce.jl:415 within &#96;_mapreduce&#39;
; │││││││┌ @ reduce.jl:257 within &#96;mapreduce_impl&#39;
          call void @j_mapreduce_impl_11840&#40;&#91;2 x double&#93;* noalias nocapture
 nonnull sret &#37;3, &#37;jl_value_t* nonnull &#37;1, i64 1, i64 &#37;9, i64 1024&#41; #0
; └└└└└└└└
  &#37;.sroa.035.0..sroa_idx39 &#61; getelementptr inbounds &#91;2 x double&#93;, &#91;2 x doub
le&#93;* &#37;3, i64 0, i64 0
  &#37;35 &#61; bitcast double* &#37;.sroa.035.0..sroa_idx39 to &lt;2 x double&gt;*
  &#37;36 &#61; load &lt;2 x double&gt;, &lt;2 x double&gt;* &#37;35, align 8
; ┌ @ reducedim.jl:719 within &#96;sum&#39;
; │┌ @ reducedim.jl:719 within &#96;#sum#627&#39;
; ││┌ @ reducedim.jl:723 within &#96;_sum&#39; @ reducedim.jl:724
; │││┌ @ reducedim.jl:310 within &#96;mapreduce&#39;
; ││││┌ @ reducedim.jl:310 within &#96;#mapreduce#620&#39;
; │││││┌ @ reducedim.jl:318 within &#96;_mapreduce_dim&#39;
; ││││││┌ @ reduce.jl:415 within &#96;_mapreduce&#39;
         br label &#37;L47

L47:                                              ; preds &#61; &#37;L43, &#37;L34, &#37;L1
7, &#37;L13
; └└└└└└└
  &#37;37 &#61; phi &lt;2 x double&gt; &#91; &#37;36, &#37;L43 &#93;, &#91; &#37;19, &#37;L13 &#93;, &#91; &#37;27, &#37;L17 &#93;, &#91; &#37;33
, &#37;L34 &#93;
; ┌ @ array.jl:219 within &#96;length&#39;
   &#37;38 &#61; bitcast &#37;jl_value_t* &#37;1 to &#37;jl_array_t*
   &#37;39 &#61; getelementptr inbounds &#37;jl_array_t, &#37;jl_array_t* &#37;38, i64 0, i32 1
   &#37;40 &#61; load i64, i64* &#37;39, align 8
; └
; ┌ @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.j
md:4 within &#96;/&#39; @ promotion.jl:314
; │┌ @ promotion.jl:282 within &#96;promote&#39;
; ││┌ @ promotion.jl:259 within &#96;_promote&#39;
; │││┌ @ number.jl:7 within &#96;convert&#39;
; ││││┌ @ float.jl:60 within &#96;Float64&#39;
       &#37;41 &#61; sitofp i64 &#37;40 to double
; │└└└└
; │ @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.j
md:4 within &#96;/&#39; @ promotion.jl:314 @ float.jl:407
   &#37;42 &#61; insertelement &lt;2 x double&gt; undef, double &#37;41, i32 0
   &#37;43 &#61; insertelement &lt;2 x double&gt; &#37;42, double &#37;41, i32 1
   &#37;44 &#61; fdiv &lt;2 x double&gt; &#37;37, &#37;43
; └
  &#37;.sroa.029.0..sroa_idx &#61; getelementptr inbounds &#91;2 x double&#93;, &#91;2 x double
&#93;* &#37;0, i64 0, i64 0
  &#37;45 &#61; bitcast double* &#37;.sroa.029.0..sroa_idx to &lt;2 x double&gt;*
  store &lt;2 x double&gt; &#37;44, &lt;2 x double&gt;* &#37;45, align 8
  ret void
&#125;
</pre>


<p>What this is doing is creating small little vectors and then parallelizing the operations of those vectors by calling specific vector-parallel instructions. Keep this in mind.</p>
<h3>Explicit SIMD</h3>
<p>The following was all a form of <strong>loop-level parallelism</strong> known as loop vectorization. It&#39;s simply easier for compilers to reason at the array level, prove iterates are independent, and automatically generate SIMD code from that. This is not necessary, and compilers can produce SIMD code from non-looping code through a process known as <strong>SLP supervectorization</strong>, but the results are far from optimal and the compiler requires a lot of time to do this calculation, meaning that it&#39;s usually not a pass used by default.</p>
<p>If you want to pack the vectors yourself, then primitives for doing so from within Julia are available in SIMD.jl. This is for &quot;real&quot; performance warriors. This looks like for example:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>SIMD</span><span class='hljl-t'>
</span><span class='hljl-n'>v</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Vec</span><span class='hljl-p'>{</span><span class='hljl-ni'>4</span><span class='hljl-p'>,</span><span class='hljl-n'>Float64</span><span class='hljl-p'>}((</span><span class='hljl-ni'>1</span><span class='hljl-p'>,</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-ni'>3</span><span class='hljl-p'>,</span><span class='hljl-ni'>4</span><span class='hljl-p'>))</span><span class='hljl-t'>
</span><span class='hljl-nd'>@show</span><span class='hljl-t'> </span><span class='hljl-n'>v</span><span class='hljl-oB'>+</span><span class='hljl-n'>v</span><span class='hljl-t'> </span><span class='hljl-cs'># basic arithmetic is supported</span><span class='hljl-t'>
</span><span class='hljl-nd'>@show</span><span class='hljl-t'> </span><span class='hljl-nf'>sum</span><span class='hljl-p'>(</span><span class='hljl-n'>v</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-cs'># basic reductions are supported</span>
</pre>


<pre class="output">
v &#43; v &#61; &lt;4 x Float64&gt;&#91;2.0, 4.0, 6.0, 8.0&#93;
sum&#40;v&#41; &#61; 10.0
10.0
</pre>


<p>Using this you can pull apart code and force the usage of SIMD vectors. One library which makes great use of this is LoopVectorization.jl. However, one word of &quot;caution&quot;:</p>
<p><strong>Most performance optimization is not trying to do something really good for performance. Most performance optimization is trying to not do something that is actively bad for performance.</strong></p>
<h3>Summary of SIMD</h3>
<ul>
<li><p>Communication in SIMD is due to locality: if things are local the processor can automatically setup the operations.</p>
</li>
<li><p>There&#39;s no real worry about &quot;getting it wrong&quot;: you cannot overwrite pieces from different parts of the arithmetic unit, and if SIMD is unsafe then it just won&#39;t auto-vectorize.</p>
</li>
<li><p>Suitable for operations measured in ns.</p>
</li>
</ul>
<h2>Next Level Up: Multithreading</h2>
<p>Last time we briefly went over multithreading and described how every process has multiple threads which share a single heap, and when multiple threads are executed simultaniously we have multithreaded parallelism. Note that you can have multiple threads which aren&#39;t executed simultaniously, like in the case of I/O operations, and this is an example of concurrency without parallelism and is commonly referred to as green threads.</p>
<p><img src="https://blog-assets.risingstack.com/2017/02/kernel-processes-and-threads-1.png" alt="" /></p>
<p>Last time we described a simple multithreaded program and noticed that multithreading has an overhead cost of around 50ns-100ns. This is due to the construction of the new stack &#40;amont other things&#41; each time a new computational thread is spun up. This means that, unlike SIMD, some thought needs to be put in as to when to perform multithreading: it&#39;s not always a good idea. It needs to be high enough on the cost for this to be counter-balanced.</p>
<p>One abstraction that was glossed over was the memory access style. Before, we were considering a single heap, or an UMA style:</p>
<p><img src="https://software.intel.com/sites/default/files/m/2/0/4/e/d/39352-figure-1.jpg" alt="" /></p>
<p>However, this is the case for all shared memory devices. For example, compute nodes on the HPC tend to be &quot;dual Xeon&quot; or &quot;quad Xeon&quot;, where each Xeon processor is itself a multicore processor. But each processor on its own accesses its own local caches, and thus one has to be aware that this is setup in a NUMA &#40;non-uniform memory access&#41; manner:</p>
<p><img src="https://software.intel.com/sites/default/files/m/2/d/c/b/2/39353-figure-2.jpg" alt="" /></p>
<p>where there is a cache that is closer to the prcessor and a cache that is further away. Care should be taken in this to localize the computation per thread, otherwise a cost associated with the memory sharing will be hit &#40;but all sharing will still be automatic&#41;.</p>
<p>In this sense, interthread communication is naturally done through the heap: if you want other threads to be able to touch a value, then you can simply place it on the heap and then it&#39;ll be available. We saw this last time by how overlapping computations can re-use the same heap-based caches, meaning that care needs to be taken with how one writes into a dynamically-allocated array.</p>
<p>A simple example that demonstrates this is. First, let&#39;s make sure we have multithreading enabled:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>Base</span><span class='hljl-oB'>.</span><span class='hljl-n'>Threads</span><span class='hljl-t'>
</span><span class='hljl-n'>Threads</span><span class='hljl-oB'>.</span><span class='hljl-nf'>nthreads</span><span class='hljl-p'>()</span><span class='hljl-t'> </span><span class='hljl-cs'># should not be 1</span>
</pre>


<pre class="output">
6
</pre>



<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>BenchmarkTools</span><span class='hljl-t'>
</span><span class='hljl-n'>acc</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-ni'>0</span><span class='hljl-t'>
</span><span class='hljl-nd'>@threads</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>10_000</span><span class='hljl-t'>
    </span><span class='hljl-kd'>global</span><span class='hljl-t'> </span><span class='hljl-n'>acc</span><span class='hljl-t'>
    </span><span class='hljl-n'>acc</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>acc</span>
</pre>


<pre class="output">
3059
</pre>


<p>The reason for this behavior is that there is a difference between the reading and the writing step to an array. Here, values are being read while other threads are writing, meaning that they see a lower value than when they are attempting to write into it. The result is that the total summation is lower than the true value because of this clashing. We can prevent this by only allowing one thread to utilize the heap-allocated variable at a time. One abstraction for doing this is <em>atomics</em>:</p>


<pre class='hljl'>
<span class='hljl-n'>acc</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Atomic</span><span class='hljl-p'>{</span><span class='hljl-n'>Int64</span><span class='hljl-p'>}(</span><span class='hljl-ni'>0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nd'>@threads</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>10_000</span><span class='hljl-t'>
    </span><span class='hljl-nf'>atomic_add!</span><span class='hljl-p'>(</span><span class='hljl-n'>acc</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>acc</span>
</pre>


<pre class="output">
Atomic&#123;Int64&#125;&#40;10000&#41;
</pre>


<p>When an atomic add is being done, all other threads wishing to do the same computation are blocked. This of course can have a massive effect on performance since atomic computations are not parallel.</p>
<p>Julia also exposes a lower level of heap control in threading using <em>locks</em></p>


<pre class='hljl'>
<span class='hljl-kd'>const</span><span class='hljl-t'> </span><span class='hljl-n'>acc_lock</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Ref</span><span class='hljl-p'>{</span><span class='hljl-n'>Int64</span><span class='hljl-p'>}(</span><span class='hljl-ni'>0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-kd'>const</span><span class='hljl-t'> </span><span class='hljl-n'>splock</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>SpinLock</span><span class='hljl-p'>()</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>f1</span><span class='hljl-p'>()</span><span class='hljl-t'>
    </span><span class='hljl-nd'>@threads</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>10_000</span><span class='hljl-t'>
        </span><span class='hljl-nf'>lock</span><span class='hljl-p'>(</span><span class='hljl-n'>splock</span><span class='hljl-p'>)</span><span class='hljl-t'>
        </span><span class='hljl-n'>acc_lock</span><span class='hljl-p'>[]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
        </span><span class='hljl-nf'>unlock</span><span class='hljl-p'>(</span><span class='hljl-n'>splock</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-kd'>const</span><span class='hljl-t'> </span><span class='hljl-n'>rsplock</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>ReentrantLock</span><span class='hljl-p'>()</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>f2</span><span class='hljl-p'>()</span><span class='hljl-t'>
    </span><span class='hljl-nd'>@threads</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>10_000</span><span class='hljl-t'>
        </span><span class='hljl-nf'>lock</span><span class='hljl-p'>(</span><span class='hljl-n'>rsplock</span><span class='hljl-p'>)</span><span class='hljl-t'>
        </span><span class='hljl-n'>acc_lock</span><span class='hljl-p'>[]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
        </span><span class='hljl-nf'>unlock</span><span class='hljl-p'>(</span><span class='hljl-n'>rsplock</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-n'>acc2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Atomic</span><span class='hljl-p'>{</span><span class='hljl-n'>Int64</span><span class='hljl-p'>}(</span><span class='hljl-ni'>0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>g</span><span class='hljl-p'>()</span><span class='hljl-t'>
  </span><span class='hljl-nd'>@threads</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>10_000</span><span class='hljl-t'>
      </span><span class='hljl-nf'>atomic_add!</span><span class='hljl-p'>(</span><span class='hljl-n'>acc2</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-p'>)</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-kd'>const</span><span class='hljl-t'> </span><span class='hljl-n'>acc_s</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Ref</span><span class='hljl-p'>{</span><span class='hljl-n'>Int64</span><span class='hljl-p'>}(</span><span class='hljl-ni'>0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>h</span><span class='hljl-p'>()</span><span class='hljl-t'>
  </span><span class='hljl-kd'>global</span><span class='hljl-t'> </span><span class='hljl-n'>acc_s</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>10_000</span><span class='hljl-t'>
      </span><span class='hljl-n'>acc_s</span><span class='hljl-p'>[]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>f1</span><span class='hljl-p'>()</span>
</pre>


<pre class="output">
460.300 μs &#40;31 allocations: 5.09 KiB&#41;
</pre>


<p><code>SpinLock</code> is non-reentrent, i.e. it will block itself if a thread that calls a <code>lock</code> does another <code>lock</code>. Therefore it has to be used with caution &#40;every <code>lock</code> goes with one <code>unlock</code>&#41;, but it&#39;s fast. <code>ReentrantLock</code> alleviates those concerns, but trades off a bit of performance:</p>


<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>f2</span><span class='hljl-p'>()</span>
</pre>


<pre class="output">
4.586 ms &#40;1984 allocations: 66.13 KiB&#41;
</pre>


<p>But if you can use atomics, they will be faster:</p>


<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>g</span><span class='hljl-p'>()</span>
</pre>


<pre class="output">
175.901 μs &#40;31 allocations: 5.09 KiB&#41;
</pre>


<p>and if your computation is actually serial, then use serial code:</p>


<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>h</span><span class='hljl-p'>()</span>
</pre>


<pre class="output">
2.344 μs &#40;0 allocations: 0 bytes&#41;
</pre>


<p>Why is this so fast? Check the code:</p>


<pre class='hljl'>
<span class='hljl-nd'>@code_llvm</span><span class='hljl-t'> </span><span class='hljl-nf'>h</span><span class='hljl-p'>()</span>
</pre>


<pre class="output">
;  @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.jm
d:26 within &#96;h&#39;
; Function Attrs: uwtable
define void @julia_h_12011&#40;&#41; #0 &#123;
top:
;  @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.jm
d:29 within &#96;h&#39;
; ┌ @ refvalue.jl:32 within &#96;getindex&#39;
; │┌ @ Base.jl:33 within &#96;getproperty&#39;
    &#37;.pre &#61; load i64, i64* inttoptr &#40;i64 931134192 to i64*&#41;, align 16
; └└
;  @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.jm
d:28 within &#96;h&#39;
  br label &#37;L2

L2:                                               ; preds &#61; &#37;L2, &#37;top
;  @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.jm
d:29 within &#96;h&#39;
; ┌ @ refvalue.jl:32 within &#96;getindex&#39;
; │┌ @ Base.jl:33 within &#96;getproperty&#39;
    &#37;0 &#61; phi i64 &#91; &#37;.pre, &#37;top &#93;, &#91; &#37;1, &#37;L2 &#93;
    &#37;value_phi &#61; phi i64 &#91; 1, &#37;top &#93;, &#91; &#37;3, &#37;L2 &#93;
; └└
; ┌ @ int.jl:86 within &#96;&#43;&#39;
   &#37;1 &#61; add i64 &#37;0, 1
; └
; ┌ @ refvalue.jl:33 within &#96;setindex&#33;&#39;
; │┌ @ Base.jl:34 within &#96;setproperty&#33;&#39;
    store i64 &#37;1, i64* inttoptr &#40;i64 931134192 to i64*&#41;, align 16
; └└
; ┌ @ range.jl:624 within &#96;iterate&#39;
; │┌ @ promotion.jl:398 within &#96;&#61;&#61;&#39;
    &#37;2 &#61; icmp eq i64 &#37;value_phi, 10000
; │└
   &#37;3 &#61; add nuw nsw i64 &#37;value_phi, 1
; └
  br i1 &#37;2, label &#37;L18, label &#37;L2

L18:                                              ; preds &#61; &#37;L2
  ret void
&#125;
</pre>


<p>It just knows to add 10000. So to get a proper timing let&#39;s make the size mutable:</p>


<pre class='hljl'>
<span class='hljl-kd'>const</span><span class='hljl-t'> </span><span class='hljl-n'>len</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Ref</span><span class='hljl-p'>{</span><span class='hljl-n'>Int</span><span class='hljl-p'>}(</span><span class='hljl-ni'>10_000</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>h2</span><span class='hljl-p'>()</span><span class='hljl-t'>
  </span><span class='hljl-kd'>global</span><span class='hljl-t'> </span><span class='hljl-n'>acc_s</span><span class='hljl-t'>
  </span><span class='hljl-kd'>global</span><span class='hljl-t'> </span><span class='hljl-n'>len</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>len</span><span class='hljl-p'>[]</span><span class='hljl-t'>
      </span><span class='hljl-n'>acc_s</span><span class='hljl-p'>[]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>h2</span><span class='hljl-p'>()</span>
</pre>


<pre class="output">
2.333 μs &#40;0 allocations: 0 bytes&#41;
</pre>



<pre class='hljl'>
<span class='hljl-nd'>@code_llvm</span><span class='hljl-t'> </span><span class='hljl-nf'>h2</span><span class='hljl-p'>()</span>
</pre>


<pre class="output">
;  @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.jm
d:3 within &#96;h2&#39;
; Function Attrs: uwtable
define void @julia_h2_12035&#40;&#41; #0 &#123;
top:
;  @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.jm
d:6 within &#96;h2&#39;
; ┌ @ refvalue.jl:32 within &#96;getindex&#39;
; │┌ @ Base.jl:33 within &#96;getproperty&#39;
    &#37;0 &#61; load i64, i64* inttoptr &#40;i64 342543248 to i64*&#41;, align 16
; └└
; ┌ @ range.jl:5 within &#96;Colon&#39;
; │┌ @ range.jl:280 within &#96;UnitRange&#39;
; ││┌ @ range.jl:285 within &#96;unitrange_last&#39;
; │││┌ @ operators.jl:350 within &#96;&gt;&#61;&#39;
; ││││┌ @ int.jl:441 within &#96;&lt;&#61;&#39;
       &#37;1 &#61; icmp sgt i64 &#37;0, 0
; │││└└
     &#37;2 &#61; select i1 &#37;1, i64 &#37;0, i64 0
; └└└
  br i1 &#37;1, label &#37;L13.preheader, label &#37;L29

L13.preheader:                                    ; preds &#61; &#37;top
;  @ C:\Users\accou\.julia\external\18337\lecture6\styles_of_parallelism.jm
d:7 within &#96;h2&#39;
; ┌ @ refvalue.jl:32 within &#96;getindex&#39;
; │┌ @ Base.jl:33 within &#96;getproperty&#39;
    &#37;.pre &#61; load i64, i64* inttoptr &#40;i64 931134192 to i64*&#41;, align 16
; └└
  br label &#37;L13

L13:                                              ; preds &#61; &#37;L13, &#37;L13.preh
eader
; ┌ @ refvalue.jl:32 within &#96;getindex&#39;
; │┌ @ Base.jl:33 within &#96;getproperty&#39;
    &#37;3 &#61; phi i64 &#91; &#37;4, &#37;L13 &#93;, &#91; &#37;.pre, &#37;L13.preheader &#93;
    &#37;value_phi2 &#61; phi i64 &#91; &#37;6, &#37;L13 &#93;, &#91; 1, &#37;L13.preheader &#93;
; └└
; ┌ @ int.jl:86 within &#96;&#43;&#39;
   &#37;4 &#61; add i64 &#37;3, 1
; └
; ┌ @ refvalue.jl:33 within &#96;setindex&#33;&#39;
; │┌ @ Base.jl:34 within &#96;setproperty&#33;&#39;
    store i64 &#37;4, i64* inttoptr &#40;i64 931134192 to i64*&#41;, align 16
; └└
; ┌ @ range.jl:624 within &#96;iterate&#39;
; │┌ @ promotion.jl:398 within &#96;&#61;&#61;&#39;
    &#37;5 &#61; icmp eq i64 &#37;value_phi2, &#37;2
; │└
   &#37;6 &#61; add nuw i64 &#37;value_phi2, 1
; └
  br i1 &#37;5, label &#37;L29, label &#37;L13

L29:                                              ; preds &#61; &#37;L13, &#37;top
  ret void
&#125;
</pre>


<p>It&#39;s still optimizing it&#33;</p>


<pre class='hljl'>
<span class='hljl-n'>non_const_len</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-ni'>10000</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>h3</span><span class='hljl-p'>()</span><span class='hljl-t'>
  </span><span class='hljl-kd'>global</span><span class='hljl-t'> </span><span class='hljl-n'>acc_s</span><span class='hljl-t'>
  </span><span class='hljl-kd'>global</span><span class='hljl-t'> </span><span class='hljl-n'>non_const_len</span><span class='hljl-t'>
  </span><span class='hljl-n'>len2</span><span class='hljl-oB'>::</span><span class='hljl-n'>Int</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>non_const_len</span><span class='hljl-t'>
  </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-n'>len2</span><span class='hljl-t'>
      </span><span class='hljl-n'>acc_s</span><span class='hljl-p'>[]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
  </span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>h3</span><span class='hljl-p'>()</span>
</pre>


<pre class="output">
2.389 μs &#40;0 allocations: 0 bytes&#41;
</pre>


<p>Note that what is shown here is a type-declaration. <code>a::T &#61; ...</code> forces <code>a</code> to be of type <code>T</code> throughout the whole function. By giving the compiler this information, I am able to use the non-constant global in a type-stable manner.</p>
<p>One last thing to note about multithreaded computations, and parallel computations, is that one cannot assume that the parallelized computation is computed in any given order. For example, the following will has a quasi-random ordering:</p>


<pre class='hljl'>
<span class='hljl-kd'>const</span><span class='hljl-t'> </span><span class='hljl-n'>a2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>zeros</span><span class='hljl-p'>(</span><span class='hljl-nf'>nthreads</span><span class='hljl-p'>()</span><span class='hljl-oB'>*</span><span class='hljl-ni'>10</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-kd'>const</span><span class='hljl-t'> </span><span class='hljl-n'>acc_lock2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Ref</span><span class='hljl-p'>{</span><span class='hljl-n'>Int64</span><span class='hljl-p'>}(</span><span class='hljl-ni'>0</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-kd'>const</span><span class='hljl-t'> </span><span class='hljl-n'>splock2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>SpinLock</span><span class='hljl-p'>()</span><span class='hljl-t'>
</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>f_order</span><span class='hljl-p'>()</span><span class='hljl-t'>
    </span><span class='hljl-nd'>@threads</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>a2</span><span class='hljl-p'>)</span><span class='hljl-t'>
        </span><span class='hljl-nf'>lock</span><span class='hljl-p'>(</span><span class='hljl-n'>splock2</span><span class='hljl-p'>)</span><span class='hljl-t'>
        </span><span class='hljl-n'>acc_lock2</span><span class='hljl-p'>[]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
        </span><span class='hljl-n'>a2</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>acc_lock2</span><span class='hljl-p'>[]</span><span class='hljl-t'>
        </span><span class='hljl-nf'>unlock</span><span class='hljl-p'>(</span><span class='hljl-n'>splock2</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>
</span><span class='hljl-nf'>f_order</span><span class='hljl-p'>()</span><span class='hljl-t'>
</span><span class='hljl-n'>a2</span>
</pre>


<pre class="output">
60-element Array&#123;Float64,1&#125;:
 41.0
 42.0
 43.0
 44.0
 45.0
 46.0
 47.0
 48.0
 49.0
 50.0
  ⋮
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0
</pre>


<p>Note that here we can see that Julia 1.5 is dividing up the work into groups of 10 for each thread, and then one thread dominates the computation at a time, but which thread dominates is random.</p>
<h3>The Dining Philosophers Problem</h3>
<p>A classic tale in parallel computing is the dining philosophers problem. In this case, there are N philosophers at a table who all want to eat at the same time, following all of the same rules. Each philosopher must alternatively think and then eat. They need both their left and right fork to start eating, but cannot start eating until they have both forks. The problem is how to setup a concurrent algorithm that will not cause any philosophers to starve.</p>
<p>The difficulty is a situation known as <em>deadlock</em>. For example, if each philosopher was told to grab the right fork when it&#39;s avaialble, and then the left fork, and put down the fork after eating, then they will all grab the right fork and none will ever eat because they will all be waiting on the left fork. This is analygous to two blocked computations which are waiting on the other to finish. Thus, when using blocking structures, one needs to be careful about deadlock&#33;</p>
<h3>Two Programming Models: Loop-Level Parallelism and Task-Based Parallelism</h3>
<p>As described in the previous lecture, one can also use <code>Threads.@spawn</code> to do multithreading in Julia v1.3&#43;. The same factors all applay: how to do locks and Mutex etc. This is a case of a parallelism construct having two alternative <strong>programming models</strong>. <code>Threads.@spawn</code> represents task-based parallelism, while <code>Threads.@threads</code> represents Loop-Level Parallelism or a parallel iterator model. Loop-based parallelization models are very high level and, assuming every iteration is independent, almost requires no code change. Task-based parallelism is a more expressive parallelism model, but usually requires modifying the code to be explicitly written as a set of parallelizable tasks. Note that in the case of Julia, <code>Threads.@threads</code> is implemented using <code>Threads.@spawn</code>&#39;s model.</p>
<h3>Summary of Multithreading</h3>
<ul>
<li><p>Communication in multithreading is done on the heap. Locks and atomics allow for a form of safe message passing.</p>
</li>
<li><p>50ns-100ns of overhead. Suitable for 1μs calculations.</p>
</li>
<li><p>Be careful of ordering and heap-allocated values.</p>
</li>
</ul>
<h2>GPU Computing</h2>
<p>GPUs are not fast. In fact, the problem with GPUs is that each processor is slow. However, GPUs have a lot of cores... like thousands.</p>
<p><img src="https://miro.medium.com/max/832/0*xzPjWMqXC0NB6D69.jpg" alt="" /></p>
<p>An RTX2080, a standard &quot;gaming&quot; GPU &#40;not even the ones in the cluster&#41;, has 2944 cores. However, not only are GPUs slow, but they also need to be programmed in a style that is <em>SPMD</em>, which standard for Single Program Multiple Data. This means that every single thread must be running the same program but on different pieces of data. Exactly the same program. If you have</p>



<pre class='hljl'>
<span class='hljl-k'>if</span><span class='hljl-t'> </span><span class='hljl-n'>a</span><span class='hljl-t'> </span><span class='hljl-oB'>&gt;</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
  </span><span class='hljl-cs'># Do something</span><span class='hljl-t'>
</span><span class='hljl-k'>else</span><span class='hljl-t'>
  </span><span class='hljl-cs'># Do something else</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span>
</pre>


<p>where some of the data goes on one branch and other data goes on the other branch, every single thread will run both branches &#40;performing &quot;fake&quot; computations while on the other branch&#41;. This means that GPU tasks should be &quot;very parallel&quot; with as few conditionals as possible.</p>
<h3>GPU Memory</h3>
<p>GPUs themselves are shared memory devices, meaning they have a heap that is shared amongst all threads. However, GPUs are heavily in the NUMA camp, where different blocks of the GPU have much faster access to certain parts of the memory. Additionally, this heap is disconnected from the standard processor, so data must be passed to the GPU and data must be returned.</p>
<p>GPU memory size is relatively small compared to CPUs. Example: the RTX2080Ti has 8GB of RAM. Thus one needs to be doing computations that are memory compact &#40;such as matrix multiplications, which are O&#40;n^3&#41; making the computation time scale quicker than the memory cost&#41;.</p>
<h3>Note on GPU Hardware</h3>
<p>Standard GPU hardware &quot;for gaming&quot;, like RTX2070, is just as fast as higher end GPU hardware for Float32. Higher end hardware, like the Tesla, add more memory, memory safety, and Float64 support. However, these require being in a server since they have alternative cooling strategies, making them a higher end product.</p>
<h3>SPMD Kernel Generation GPU Computing Models</h3>
<p>The core programming models for GPU computing are SPMD kernel compilers, of which the most well-known is CUDA. CUDA is a C&#43;&#43;-like programming language which compiles to .ptx kernels, and GPU execution on NVIDIA GPUs is done by &quot;all steams&quot; of a GPU doing concurrent execution of the kernel &#40;generally, without going into more details, you can of &quot;all streams&quot; as just meaning &quot;all cores&quot;. More detailed views of GPU execution will come later&#41;.</p>
<p>.ptx CUDA kernels can be compiled from LLVM IR, and thus since Julia is a programming language which emits LLVM IR for all of its operations, native Julia programs are compatible with compilation to CUDA. The helper functions to enable this separate compilation path is CUDA.jl. Let&#39;s take a look at a basic CUDA.jl kernel generating example:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>CUDA</span><span class='hljl-t'>

</span><span class='hljl-n'>N</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-ni'>2</span><span class='hljl-oB'>^</span><span class='hljl-ni'>20</span><span class='hljl-t'>
</span><span class='hljl-n'>x_d</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>CUDA</span><span class='hljl-oB'>.</span><span class='hljl-nf'>fill</span><span class='hljl-p'>(</span><span class='hljl-nfB'>1.0f0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>N</span><span class='hljl-p'>)</span><span class='hljl-t'>  </span><span class='hljl-cs'># a vector stored on the GPU filled with 1.0 (Float32)</span><span class='hljl-t'>
</span><span class='hljl-n'>y_d</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>CUDA</span><span class='hljl-oB'>.</span><span class='hljl-nf'>fill</span><span class='hljl-p'>(</span><span class='hljl-nfB'>2.0f0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>N</span><span class='hljl-p'>)</span><span class='hljl-t'>  </span><span class='hljl-cs'># a vector stored on the GPU filled with 2.0</span><span class='hljl-t'>

</span><span class='hljl-k'>function</span><span class='hljl-t'> </span><span class='hljl-nf'>gpu_add2!</span><span class='hljl-p'>(</span><span class='hljl-n'>y</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>x</span><span class='hljl-p'>)</span><span class='hljl-t'>
    </span><span class='hljl-n'>index</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>threadIdx</span><span class='hljl-p'>()</span><span class='hljl-oB'>.</span><span class='hljl-n'>x</span><span class='hljl-t'>    </span><span class='hljl-cs'># this example only requires linear indexing, so just use `x`</span><span class='hljl-t'>
    </span><span class='hljl-n'>stride</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>blockDim</span><span class='hljl-p'>()</span><span class='hljl-oB'>.</span><span class='hljl-n'>x</span><span class='hljl-t'>
    </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>index</span><span class='hljl-oB'>:</span><span class='hljl-n'>stride</span><span class='hljl-oB'>:</span><span class='hljl-nf'>length</span><span class='hljl-p'>(</span><span class='hljl-n'>y</span><span class='hljl-p'>)</span><span class='hljl-t'>
        </span><span class='hljl-nd'>@inbounds</span><span class='hljl-t'> </span><span class='hljl-n'>y</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>]</span><span class='hljl-t'> </span><span class='hljl-oB'>+=</span><span class='hljl-t'> </span><span class='hljl-n'>x</span><span class='hljl-p'>[</span><span class='hljl-n'>i</span><span class='hljl-p'>]</span><span class='hljl-t'>
    </span><span class='hljl-k'>end</span><span class='hljl-t'>
    </span><span class='hljl-k'>return</span><span class='hljl-t'> </span><span class='hljl-n'>nothing</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span><span class='hljl-t'>

</span><span class='hljl-nf'>fill!</span><span class='hljl-p'>(</span><span class='hljl-n'>y_d</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-ni'>2</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nd'>@cuda</span><span class='hljl-t'> </span><span class='hljl-n'>threads</span><span class='hljl-oB'>=</span><span class='hljl-ni'>256</span><span class='hljl-t'> </span><span class='hljl-nf'>gpu_add2!</span><span class='hljl-p'>(</span><span class='hljl-n'>y_d</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>x_d</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>all</span><span class='hljl-p'>(</span><span class='hljl-nf'>Array</span><span class='hljl-p'>(</span><span class='hljl-n'>y_d</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>.==</span><span class='hljl-t'> </span><span class='hljl-nfB'>3.0f0</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
true
</pre>


<p>The key to understanding the SPMD kernel approach is the <code>index &#61; threadIdx&#40;&#41;.x</code> and <code>stride &#61; blockDim&#40;&#41;.x</code> portions.</p>
<p><img src="https://juliagpu.gitlab.io/CUDA.jl/tutorials/intro1.png" alt="" /></p>
<p>The way kernels are expected to run in parallel is that they are given a specific block of the computation and are expected to write a kernel which only on that small block of the input. This kernel is then called on every separate thread on the GPU, making each CUDA core simultaniously compute each block. Thus as a user in such a SPMD programming model, you never specify the computation globally but instead simply specify how chunks should behave, giving the compiler the leeway to determine the optimal global execution.</p>
<h3>Array-Based GPU Computing Models</h3>
<p>The simplest version of GPU computing is the array-based programming model.</p>


<pre class='hljl'>
<span class='hljl-n'>A</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>rand</span><span class='hljl-p'>(</span><span class='hljl-ni'>100</span><span class='hljl-p'>,</span><span class='hljl-ni'>100</span><span class='hljl-p'>);</span><span class='hljl-t'> </span><span class='hljl-n'>B</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>rand</span><span class='hljl-p'>(</span><span class='hljl-ni'>100</span><span class='hljl-p'>,</span><span class='hljl-ni'>100</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>CUDA</span><span class='hljl-t'>
</span><span class='hljl-cs'># Pass to the GPU</span><span class='hljl-t'>
</span><span class='hljl-n'>cuA</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>cu</span><span class='hljl-p'>(</span><span class='hljl-n'>A</span><span class='hljl-p'>);</span><span class='hljl-t'> </span><span class='hljl-n'>cuB</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>cu</span><span class='hljl-p'>(</span><span class='hljl-n'>B</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>cuC</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>cuA</span><span class='hljl-oB'>*</span><span class='hljl-n'>cuB</span><span class='hljl-t'>
</span><span class='hljl-cs'># Pass to the CPU</span><span class='hljl-t'>
</span><span class='hljl-n'>C</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Array</span><span class='hljl-p'>(</span><span class='hljl-n'>cuC</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
100×100 Array&#123;Float32,2&#125;:
 21.3225  25.4278  24.5773  26.6089  …  21.0833  25.0389  22.7444  25.5705
 25.8577  27.2815  28.2439  31.2978     24.4566  27.6241  25.8184  30.2117
 23.2739  23.9996  26.5368  27.0796     22.6042  23.3092  23.689   25.9774
 23.6532  25.7988  25.0622  27.7568     21.9137  26.1695  24.0751  29.1244
 21.8137  25.4429  24.6081  27.8353     22.1399  24.7156  22.8079  25.912
 22.8847  26.2641  25.951   27.4366  …  22.5984  26.7368  23.2251  25.9416
 23.2013  23.8571  24.165   27.7386     20.6205  24.987   22.2018  26.1901
 25.7156  26.6828  28.2383  30.6743     23.2989  29.3801  25.5503  29.0319
 23.7644  22.2289  24.0214  25.6603     21.4567  25.3679  23.1681  26.6249
 24.1221  23.8238  25.4609  27.4193     21.5165  24.0425  23.123   28.4186
  ⋮                                  ⋱                             
 25.0148  26.0299  27.0189  28.7919     23.7598  26.5612  24.9481  27.4216
 24.1719  25.0543  25.8945  25.3179     21.8304  24.8503  22.2669  27.7159
 24.3072  24.9768  27.0445  29.6013     23.1415  26.5856  25.643   27.3642
 23.6446  26.3437  27.365   28.5331     23.3952  27.2082  24.6372  26.6592
 23.6974  23.2991  24.4604  26.6876  …  20.4468  22.3201  22.0291  26.4612
 25.6882  27.7116  27.6376  31.5592     23.8116  27.0794  25.8074  27.6167
 21.9317  24.277   26.2441  27.3912     21.1978  23.5688  22.7829  26.5278
 27.9257  26.9172  27.7042  30.7314     23.5327  29.8594  27.7334  31.6945
 24.9751  27.5279  28.7485  29.3244     24.2676  27.5577  25.1622  28.8218
</pre>


<p>Let&#39;s see the transfer times:</p>


<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>cu</span><span class='hljl-p'>(</span><span class='hljl-n'>A</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
28.100 μs &#40;10 allocations: 39.33 KiB&#41;
100×100 CuArray&#123;Float32,2&#125;:
 0.139109   0.782002   0.457874   …  0.119527  0.286168   0.509089
 0.999584   0.17679    0.168456      0.916752  0.315931   0.596657
 0.734683   0.0984683  0.049453      0.65644   0.488426   0.236935
 0.678611   0.653062   0.180714      0.541439  0.628291   0.0742101
 0.411425   0.313711   0.0446392     0.897123  0.41332    0.239645
 0.483069   0.549755   0.430024   …  0.521924  0.827396   0.591979
 0.0922345  0.803593   0.41104       0.377889  0.540939   0.32528
 0.304025   0.692909   0.895567      0.361009  0.715111   0.586864
 0.807941   0.873825   0.936103      0.407117  0.67271    0.0193138
 0.498051   0.789486   0.330725      0.599514  0.301894   0.371843
 ⋮                                ⋱                       
 0.97266    0.870757   0.402816      0.832312  0.942345   0.214438
 0.0232706  0.506793   0.469258      0.728999  0.786748   0.111056
 0.254682   0.761862   0.844721      0.566113  0.0579659  0.937348
 0.692662   0.232629   0.486957      0.809011  0.707646   0.143993
 0.0466165  0.319511   0.280374   …  0.82787   0.80973    0.670359
 0.140467   0.744883   0.793176      0.561578  0.0205778  0.762718
 0.364728   0.291075   0.390448      0.437533  0.360211   0.166104
 0.0120339  0.741594   0.564586      0.786979  0.268043   0.337639
 0.527659   0.912363   0.447123      0.881737  0.410785   0.31833
</pre>



<pre class='hljl'>
<span class='hljl-nd'>@btime</span><span class='hljl-t'> </span><span class='hljl-nf'>Array</span><span class='hljl-p'>(</span><span class='hljl-n'>cuC</span><span class='hljl-p'>)</span>
</pre>


<pre class="output">
44.400 μs &#40;4 allocations: 39.17 KiB&#41;
100×100 Array&#123;Float32,2&#125;:
 21.3225  25.4278  24.5773  26.6089  …  21.0833  25.0389  22.7444  25.5705
 25.8577  27.2815  28.2439  31.2978     24.4566  27.6241  25.8184  30.2117
 23.2739  23.9996  26.5368  27.0796     22.6042  23.3092  23.689   25.9774
 23.6532  25.7988  25.0622  27.7568     21.9137  26.1695  24.0751  29.1244
 21.8137  25.4429  24.6081  27.8353     22.1399  24.7156  22.8079  25.912
 22.8847  26.2641  25.951   27.4366  …  22.5984  26.7368  23.2251  25.9416
 23.2013  23.8571  24.165   27.7386     20.6205  24.987   22.2018  26.1901
 25.7156  26.6828  28.2383  30.6743     23.2989  29.3801  25.5503  29.0319
 23.7644  22.2289  24.0214  25.6603     21.4567  25.3679  23.1681  26.6249
 24.1221  23.8238  25.4609  27.4193     21.5165  24.0425  23.123   28.4186
  ⋮                                  ⋱                             
 25.0148  26.0299  27.0189  28.7919     23.7598  26.5612  24.9481  27.4216
 24.1719  25.0543  25.8945  25.3179     21.8304  24.8503  22.2669  27.7159
 24.3072  24.9768  27.0445  29.6013     23.1415  26.5856  25.643   27.3642
 23.6446  26.3437  27.365   28.5331     23.3952  27.2082  24.6372  26.6592
 23.6974  23.2991  24.4604  26.6876  …  20.4468  22.3201  22.0291  26.4612
 25.6882  27.7116  27.6376  31.5592     23.8116  27.0794  25.8074  27.6167
 21.9317  24.277   26.2441  27.3912     21.1978  23.5688  22.7829  26.5278
 27.9257  26.9172  27.7042  30.7314     23.5327  29.8594  27.7334  31.6945
 24.9751  27.5279  28.7485  29.3244     24.2676  27.5577  25.1622  28.8218
</pre>


<p>The cost transferring is about 20μs-50μs in each direction, meaning that one needs to be doing operations that cost at least 200μs for GPUs to break even. A good rule of thumb is that GPU computations should take at least a milisecond, or GPU memory should be re-used.</p>
<h3>Summary of GPUs</h3>
<ul>
<li><p>GPUs cores are slow</p>
</li>
<li><p>GPUs are SPMD</p>
</li>
<li><p>GPUs are generally used for linear algebra</p>
</li>
<li><p>Suitable for SPMD 1ms computations</p>
</li>
</ul>
<h2>Xeon Phi Accelerators and OpenCL</h2>
<p>Other architectures exist to keep in mind. Xeon Phis are a now-defunct accelerator that used X86 &#40;standard processors&#41; as the base, using hundreds of them. For example, the Knights Landing series had 256 core accelerator cards. These were all clocked down, meaning they were still slower than a standard CPU, but there were less restrictions on SPMD &#40;though SPMD-like computations were still preferred in order to heavily make use of SIMD&#41;. However, because machine learning essentially only needs linear algebra, and linear algebra is faster when restricting to SPMD-architectures, this failed. These devices can still be found on many high end clusters.</p>
<p>One alternative to CUDA is OpenCL which supports alternative architectures such as the Xeon Phi at the same time that it supports GPUs. However, one of the issues with OpenCL is that its BLAS implementation currently does not match the speed of CuBLAS, which makes NVIDIA-specific libraries still the king of machine learning and most scientific computing.</p>
<h2>TPU Computing</h2>
<p>TPUs are tensor processing units, which is Google&#39;s newest accelerator technology. They are essentially just &quot;tensor operation compilers&quot;, which in computer science speak is simply higher dimensional linear algebra. To do this, they internally utilize a BFloat16 type, which is a 16-bit floating point number with the same exponent size as a Float32 with an 8-bit significand. This means that computations are highly prone to <em>catastrophic cancellation</em>. This computational device only works because BFloat16 has primitive operations for FMA which allows 32-bit-like accuracy of multiply-add operations, and thus computations which are only dot products &#40;linear algebra&#41; end up okay. Thus this is simply a GPU-like device which has gone further to completely specialize in linear algebra.</p>
<h2>Multiprocessing &#40;Distributed Computing&#41;</h2>
<p>While multithreading computes with multiple threads, multiprocessing computes with multiple independent processes. Note that processes do not share any memory, not heap or data, and thus this mode of computing also allows for <em>distributed computations</em>, which is the case where processes may be on separate computing hardware. However, even if they are on the same hardware, the lack of a shared address space means that multiprocessing has to do <em>message passing</em>, i.e. send data from one process to the other.</p>
<h3>Distributed Tasks with Explicit Memory Handling: The Master-Worker Model</h3>
<p>Given the amount of control over data handling, there are many different models for distributed computing. The simplest, the one that Julia&#39;s Distributed Standard Library defaults to, is the <em>master-worker model</em>. The master-worker model has one process, deemed the master, which controls the worker processes.</p>
<p>Here we can start by adding some new worker processes:</p>



<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>Distributed</span><span class='hljl-t'>
</span><span class='hljl-nf'>addprocs</span><span class='hljl-p'>(</span><span class='hljl-ni'>4</span><span class='hljl-p'>)</span>
</pre>


<p>This adds 4 worker processes for the master to control. The simplest computations are those where the master process gives the worker process a job which returns the value afterwards. For example, a <code>pmap</code> operation or <code>@distributed</code> loop gives the worker a function to execute, along with the data, and the worker then computes and returns the result.</p>
<p>At a lower level, this is done by <code>Distributed.@spawn</code>ing jobs, or using a <code>remotecall</code> and <code>fetch</code>ing the result. <a href="https://github.com/ChrisRackauckas/ParallelDataTransfer.jl">ParallelDataTransfer.jl</a> gives an extended set of primitive message passing operations. For example, we can explicity tell it to compute a function <code>f</code> on the remote process like:</p>



<pre class='hljl'>
<span class='hljl-nd'>@everywhere</span><span class='hljl-t'> </span><span class='hljl-nf'>f</span><span class='hljl-p'>(</span><span class='hljl-n'>x</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>x</span><span class='hljl-oB'>.^</span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-cs'># Define this function on all processes</span><span class='hljl-t'>
</span><span class='hljl-n'>t</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>remotecall</span><span class='hljl-p'>(</span><span class='hljl-n'>f</span><span class='hljl-p'>,</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-nf'>randn</span><span class='hljl-p'>(</span><span class='hljl-ni'>10</span><span class='hljl-p'>))</span>
</pre>


<p><code>remotecall</code> is a non-blocking operation that returns a <code>Future</code>. To access the data, one should use the blocking operation <code>fetch</code> to receive the data:</p>



<pre class='hljl'>
<span class='hljl-n'>xsq</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>fetch</span><span class='hljl-p'>(</span><span class='hljl-n'>t</span><span class='hljl-p'>)</span>
</pre>


<h3>Distributed Tasks with Implicit Memory Handling: Distributed Task-Based Parallelism</h3>
<p>Another popular programming model for distributed computation is task-based parallelism but where all of the memory handling is implicit. Since, unlike the shared memory parallelism case, data transfers are required for given processes to share heap allocated values, distributed task-based parallelism libraries tend to want a global view of the whole computation in order to build a sophisticated schedule that includes where certain data lives and when transfers will occur. Because of this, distributed task-based parallelism libraries tend to want the entire <strong>computational graph</strong> of the computation, to be able to restructure the graph as necessary with their own data transfer portions spliced into the compute. Examples of this kind of framework are:</p>
<ul>
<li><p>Tensorflow</p>
</li>
<li><p>dask &#40;&quot;distributed tasks&quot;&#41;</p>
</li>
<li><p>Dagger.jl</p>
</li>
</ul>
<p>Using these kinds of libraries requires building a directed acyclic graph &#40;DAG&#41;. For example, the following showcases how to use Dagger.jl to reprsent a bunch of summations:</p>



<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>Dagger</span><span class='hljl-t'>

</span><span class='hljl-nf'>add1</span><span class='hljl-p'>(</span><span class='hljl-n'>value</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>value</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'>
</span><span class='hljl-nf'>add2</span><span class='hljl-p'>(</span><span class='hljl-n'>value</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>value</span><span class='hljl-t'> </span><span class='hljl-oB'>+</span><span class='hljl-t'> </span><span class='hljl-ni'>2</span><span class='hljl-t'>
</span><span class='hljl-nf'>combine</span><span class='hljl-p'>(</span><span class='hljl-n'>a</span><span class='hljl-oB'>...</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>sum</span><span class='hljl-p'>(</span><span class='hljl-n'>a</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>delayed</span><span class='hljl-p'>(</span><span class='hljl-n'>add1</span><span class='hljl-p'>)(</span><span class='hljl-ni'>4</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>q</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>delayed</span><span class='hljl-p'>(</span><span class='hljl-n'>add2</span><span class='hljl-p'>)(</span><span class='hljl-n'>p</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>r</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>delayed</span><span class='hljl-p'>(</span><span class='hljl-n'>add1</span><span class='hljl-p'>)(</span><span class='hljl-ni'>3</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>s</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>delayed</span><span class='hljl-p'>(</span><span class='hljl-n'>combine</span><span class='hljl-p'>)(</span><span class='hljl-n'>p</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>q</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>r</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-nd'>@assert</span><span class='hljl-t'> </span><span class='hljl-nf'>collect</span><span class='hljl-p'>(</span><span class='hljl-n'>s</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-oB'>==</span><span class='hljl-t'> </span><span class='hljl-ni'>16</span>
</pre>


<p>Once the global computation is specified, commands like <code>collect</code> are used to instantiate the graph on given input data, which then run the computation in a &#40;potentially&#41; distributed manner, depending on internal scheduler heuristics.</p>
<h3>Distributed Array-Based Parallelism: SharedArrays, Elemental, and DArrays</h3>
<p>Because array operations are a standard way to compute in scientific computing, there are higher level primitives to help with message passing. A <code>SharedArray</code> is an array which acts like a shared memory device. This means that every change to a <code>SharedArray</code> causes message passing to keep them in sync, and thus this should be used with a performance caution. <a href="https://github.com/JuliaParallel/DistributedArrays.jl">DistributedArrays.jl</a> is a parallel array type which has local blocks and can be used for writing higher level abstractions with explicit message passing. Because it is currently missing high-level parallel linear algebra, currently the recommended tool for distributed linear algebra is <a href="https://github.com/JuliaParallel/Elemental.jl">Elemental.jl</a>.</p>
<h3>MapReduce, Hadoop, and Spark: The Map-Reduce Model</h3>
<p>Many data-parallel operations work by mapping a function <code>f</code> onto each piece of data and then reducing it. For example, the sum of squares maps the function <code>x -&gt; x^2</code> onto each value, and then these values are reduced by performing a summation. MapReduce was a Google framework in the 2000&#39;s built around this as the parallel computing concept, and current data-handling frameworks, like Hadoop and Spark, continue this as the core distributed programming model.</p>
<p>In Julia, there exists the <code>mapreduce</code> function for performing serial mapreduce operations. It also work on GPUs. However, it does not auto-distribute. For distributed map-reduce programming, the <code>@distributed</code> for-loop macro can be used. For example, sum of squares of random numbers is:</p>



<pre class='hljl'>
<span class='hljl-nd'>@distributed</span><span class='hljl-t'> </span><span class='hljl-p'>(</span><span class='hljl-oB'>+</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-kp'>in</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>1000</span><span class='hljl-t'>
  </span><span class='hljl-nf'>rand</span><span class='hljl-p'>()</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-t'>
</span><span class='hljl-k'>end</span>
</pre>


<p>One can see that computing summary statistics is easily done in this framework which is why it was majorly adopted among &quot;big data&quot; communities.</p>
<p><code>@distributed</code> uses a static scheduler. The dynamic scheduling equivalent is <code>pmap</code>:</p>



<pre class='hljl'>
<span class='hljl-nf'>pmap</span><span class='hljl-p'>(</span><span class='hljl-n'>i</span><span class='hljl-oB'>-&gt;</span><span class='hljl-nf'>rand</span><span class='hljl-p'>()</span><span class='hljl-oB'>^</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>100</span><span class='hljl-p'>)</span>
</pre>


<p>which will dynamically allocate jobs to processes as they declare they have finished jobs. This thus has the same performance difference behavior as <code>Threads.@threads</code> vs <code>Threads.@spawn</code>.</p>
<h3>MPI: The Distributed SPMD Model</h3>
<p>The main way to do high-performace multiprocessing is <em>MPI</em>, which is an old distributed computing interface from the C/Fortran days. Julia has access to the MPI programming model through MPI.jl. The programming model for MPI is that every computer is running the same program, and synchronization is performed by blocking communication. For example, let&#39;s look at the following:</p>



<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>MPI</span><span class='hljl-t'>
</span><span class='hljl-n'>MPI</span><span class='hljl-oB'>.</span><span class='hljl-nf'>Init</span><span class='hljl-p'>()</span><span class='hljl-t'>

</span><span class='hljl-n'>comm</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>MPI</span><span class='hljl-oB'>.</span><span class='hljl-n'>COMM_WORLD</span><span class='hljl-t'>
</span><span class='hljl-n'>rank</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>MPI</span><span class='hljl-oB'>.</span><span class='hljl-nf'>Comm_rank</span><span class='hljl-p'>(</span><span class='hljl-n'>comm</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>size</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>MPI</span><span class='hljl-oB'>.</span><span class='hljl-nf'>Comm_size</span><span class='hljl-p'>(</span><span class='hljl-n'>comm</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-n'>dst</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>mod</span><span class='hljl-p'>(</span><span class='hljl-n'>rank</span><span class='hljl-oB'>+</span><span class='hljl-ni'>1</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>size</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>src</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>mod</span><span class='hljl-p'>(</span><span class='hljl-n'>rank</span><span class='hljl-oB'>-</span><span class='hljl-ni'>1</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>size</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-n'>N</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-ni'>4</span><span class='hljl-t'>

</span><span class='hljl-n'>send_mesg</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Array</span><span class='hljl-p'>{</span><span class='hljl-n'>Float64</span><span class='hljl-p'>}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>N</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>recv_mesg</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>Array</span><span class='hljl-p'>{</span><span class='hljl-n'>Float64</span><span class='hljl-p'>}(</span><span class='hljl-n'>undef</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>N</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-nf'>fill!</span><span class='hljl-p'>(</span><span class='hljl-n'>send_mesg</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nf'>Float64</span><span class='hljl-p'>(</span><span class='hljl-n'>rank</span><span class='hljl-p'>))</span><span class='hljl-t'>

</span><span class='hljl-n'>rreq</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>MPI</span><span class='hljl-oB'>.</span><span class='hljl-nf'>Irecv!</span><span class='hljl-p'>(</span><span class='hljl-n'>recv_mesg</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>src</span><span class='hljl-p'>,</span><span class='hljl-t'>  </span><span class='hljl-n'>src</span><span class='hljl-oB'>+</span><span class='hljl-ni'>32</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>comm</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-nf'>print</span><span class='hljl-p'>(</span><span class='hljl-s'>&quot;</span><span class='hljl-si'>$rank</span><span class='hljl-s'>: Sending   </span><span class='hljl-si'>$rank</span><span class='hljl-s'> -&gt; </span><span class='hljl-si'>$dst</span><span class='hljl-s'> = </span><span class='hljl-si'>$send_mesg</span><span class='hljl-se'>\n</span><span class='hljl-s'>&quot;</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>sreq</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>MPI</span><span class='hljl-oB'>.</span><span class='hljl-nf'>Isend</span><span class='hljl-p'>(</span><span class='hljl-n'>send_mesg</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>dst</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>rank</span><span class='hljl-oB'>+</span><span class='hljl-ni'>32</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>comm</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-n'>stats</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-n'>MPI</span><span class='hljl-oB'>.</span><span class='hljl-nf'>Waitall!</span><span class='hljl-p'>([</span><span class='hljl-n'>rreq</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>sreq</span><span class='hljl-p'>])</span><span class='hljl-t'>

</span><span class='hljl-nf'>print</span><span class='hljl-p'>(</span><span class='hljl-s'>&quot;</span><span class='hljl-si'>$rank</span><span class='hljl-s'>: Received </span><span class='hljl-si'>$src</span><span class='hljl-s'> -&gt; </span><span class='hljl-si'>$rank</span><span class='hljl-s'> = </span><span class='hljl-si'>$recv_mesg</span><span class='hljl-se'>\n</span><span class='hljl-s'>&quot;</span><span class='hljl-p'>)</span><span class='hljl-t'>

</span><span class='hljl-n'>MPI</span><span class='hljl-oB'>.</span><span class='hljl-nf'>Barrier</span><span class='hljl-p'>(</span><span class='hljl-n'>comm</span><span class='hljl-p'>)</span>
</pre>




<pre class='hljl'>
<span class='hljl-oB'>&gt;</span><span class='hljl-t'> </span><span class='hljl-n'>mpiexecjl</span><span class='hljl-t'> </span><span class='hljl-oB'>-</span><span class='hljl-n'>n</span><span class='hljl-t'> </span><span class='hljl-ni'>3</span><span class='hljl-t'> </span><span class='hljl-n'>julia</span><span class='hljl-t'> </span><span class='hljl-n'>examples</span><span class='hljl-oB'>/</span><span class='hljl-ni'>04</span><span class='hljl-oB'>-</span><span class='hljl-n'>sendrecv</span><span class='hljl-oB'>.</span><span class='hljl-n'>jl</span><span class='hljl-t'>
</span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-t'> </span><span class='hljl-n'>Sending</span><span class='hljl-t'>   </span><span class='hljl-ni'>1</span><span class='hljl-t'> </span><span class='hljl-oB'>-&gt;</span><span class='hljl-t'> </span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-ni'>0</span><span class='hljl-oB'>:</span><span class='hljl-t'> </span><span class='hljl-n'>Sending</span><span class='hljl-t'>   </span><span class='hljl-ni'>0</span><span class='hljl-t'> </span><span class='hljl-oB'>-&gt;</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-t'> </span><span class='hljl-n'>Received</span><span class='hljl-t'> </span><span class='hljl-ni'>0</span><span class='hljl-t'> </span><span class='hljl-oB'>-&gt;</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>0.0</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-ni'>2</span><span class='hljl-oB'>:</span><span class='hljl-t'> </span><span class='hljl-n'>Sending</span><span class='hljl-t'>   </span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-oB'>-&gt;</span><span class='hljl-t'> </span><span class='hljl-ni'>0</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-ni'>0</span><span class='hljl-oB'>:</span><span class='hljl-t'> </span><span class='hljl-n'>Received</span><span class='hljl-t'> </span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-oB'>-&gt;</span><span class='hljl-t'> </span><span class='hljl-ni'>0</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>2.0</span><span class='hljl-p'>]</span><span class='hljl-t'>
</span><span class='hljl-ni'>2</span><span class='hljl-oB'>:</span><span class='hljl-t'> </span><span class='hljl-n'>Received</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-t'> </span><span class='hljl-oB'>-&gt;</span><span class='hljl-t'> </span><span class='hljl-ni'>2</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-p'>[</span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-nfB'>1.0</span><span class='hljl-p'>]</span>
</pre>


<p>Let&#39;s investigate this a little bit. Think about having two computers run this line-by-line side by side. They will both locally build arrays, and then call <code>MPI.Irecv&#33;</code>, which is an asynchronous non-blocking call to listen for a message from a given <code>rank</code> &#40;a rank is the ID for a given process&#41;. Then they call their <code>sreq &#61; MPI.Isend</code> function, which is an asynchronous non-blocking call to send a message <code>send_mesg</code> to the chosen <code>rank</code>. When the expected message is found, <code>MPI.Irecv&#33;</code> will then run on its green thread and finish, updating the <code>recv_mesg</code> with the information from the message. However, in order to make sure all of the messages are received, we have added in a blocking operation <code>MPI.Waitall&#33;&#40;&#91;rreq, sreq&#93;&#41;</code>, which will block all further execution on the given rank until both its <code>rreq</code> and <code>sreq</code> tasks are completed. After that is done, each given rank will have its updated data, and the script will continue on all ranks.</p>
<p>This model is thus very asynchronous and allows for many different computers to run one highly parallelized program, managing the data transmissions in a sparse way without a single computer in charge of managing the whole computation. However, it can be prone to deadlock, since errors in the program may for example require rank 1 to receive a message from rank 2 before continuing the program, but rank 2 won&#39;t continue to program until it receives a message from rank 1. For this reason, while MPI has been the most successful large-scale distributed computing model and almost all major high-performance computing &#40;HPC&#41; cluster competitions have been won by codes utilizing the MPI model, the MPI model is nowadays considered a last resort due to these safety issues.</p>
<h3>Summary of Multiprocessing</h3>
<ul>
<li><p>Cost is hardware dependent: only suitable for 1ms or higher depending on the connections through which the messages are being passed and the topology of the network.</p>
</li>
<li><p>The Master-worker programming model is Julia&#39;s <code>Distributed</code> model</p>
</li>
<li><p>The Map-reduce programming model is a common data-handling model</p>
</li>
<li><p>Array-based distributed computations are another abstraction, used in all forms of parallelism.</p>
</li>
<li><p>MPI is a SPMD model of distributed computing, where each process is completely independent and one just controls the memory handling.</p>
</li>
</ul>
<h2>The Bait-and-switch: Parallelism is about Programming Models</h2>
<p>While this looked like a lecture about parallel programming at the different levels and types of hardware, this wide overview showcases that the real underlying commonality within parallel program is in the <strong>parallel programming models</strong>, of which there are not too many. There are:</p>
<ul>
<li><p>Map-reduce parallelism models. <code>pmap</code>, MapReduce &#40;Hadoop/Spark&#41;</p>
<ul>
<li><p>Pros: Easy to use</p>
</li>
<li><p>Cons: Requires that your program is specifically only mapping functions <code>f</code> and reducing them. That said, many data science operations like <code>mean</code>, <code>variance</code>, <code>maximum</code>, etc. can be represented as map-reduce calls, which lead to the popularity of these approaches for &quot;big data&quot; operations.</p>
</li>
</ul>
</li>
<li><p>Array-based parallelism models. SIMD &#40;at the compiler level&#41;, <code>CuArray</code>, <code>DistributedArray</code>, <code>PyTorch.torch</code>, ...</p>
<ul>
<li><p>Pros: Easy to use, can have very fast library implementations for specific functions</p>
</li>
<li><p>Cons: Less control and restricted to specific functions implemented by the library. Parallelism matches the data structure, so it requires the user to be careful and know the best way to split the data.</p>
</li>
</ul>
</li>
<li><p>Loop-based parallelism models. <code>Threads.@threads</code>, <code>@distributed</code>, OpenMP, MATLAB&#39;s <code>parfor</code>, Chapel&#39;s iterator parallelism</p>
<ul>
<li><p>Pros: Easy to use, almost no code change can make existing loops parallelized</p>
</li>
<li><p>Cons: Refined operations, like locking and sharing data, can be awkward to write. Less control over fine details like scheduling, meaning less opportunities to optimize.</p>
</li>
</ul>
</li>
<li><p>Task-based parallelism models with implicit distributed data handling. <code>Threads.@spawn</code>, Dagger.jl, TensorFlow, dask</p>
<ul>
<li><p>Pros: Relatively high level, low risk of errors since parallelism is mostly handled for the user. User simply describes which functions to call in what order.</p>
</li>
<li><p>Cons: When used on distributed systems, implicit data handling is hard, meaning it&#39;s generally not as efficient if you don&#39;t optimize the code yourself or help the optimizer, and these require specific programming constructs for building the computational graph. Note this is only a downside for distributed data parallelism, whereas when applied to shared memory systems these aspects no longer require handling by the task scheduler.</p>
</li>
</ul>
</li>
<li><p>Task-based parallelism models with explicit data handling. <code>Distributed.@spawn</code></p>
<ul>
<li><p>Pros: Allows for control over what compute hardware will have specific pieces of data and allows for transfering data manually.</p>
</li>
<li><p>Cons: Requires transfering data manually. All computations are managed by a single process/computer/node and thus it can have some issues scaling to extreme &#40;1000&#43; node&#41; computing situtations.</p>
</li>
</ul>
</li>
<li><p>SPMD kernel parallelism models. CUDA, MPI, KernelAbstractions.jl</p>
<ul>
<li><p>Pros: Reduces the problem for the user to only specify what happens in small chunks of the problem. Works on accelerator hardware like GPUs, TPUs, and beyond.</p>
</li>
<li><p>Cons: Only works for computations that be represented block-wise, and relies on the compiler to generate good code.</p>
</li>
</ul>
</li>
</ul>
<p>In this sense, the different parallel programming &quot;languages&quot; and features are much more similar than they are all different, falling into similar categories.</p>


        <HR/>
        <div class="footer">
          <p>
            Published from <a href="styles_of_parallelism.jmd">styles_of_parallelism.jmd</a>
            using <a href="http://github.com/JunoLab/Weave.jl">Weave.jl</a> v0.10.3 on 2020-09-25.
          </p>
        </div>
      </div>
    </div>
  </div>
</BODY>

</HTML>
